clear; clc
%% create model
g85=2.0*pi* 4.7; g87=2.0*pi* 7.0; % 2pi*Hz/nT

R0 = 10.0e3; % s^{-1}
T1 = 1.0e-3; T2 = 15.0e-6; %s

fc = 34.0e3; % Hz

bx=20.0; by = 0.0; %nT 

ratio = 0.5;
mod_amp = ratio * 2*pi*fc; mod_offset = 2*pi*fc;


m = Algorithm.BlochEquation(T1, T2, 'solution', 'periodic', 'period', 1/fc, 'steady', [0.0 0.0 1.0], 'steadyOpt', 'manual');


m.set_Omega('omegaZ', m.cosine_pump(mod_amp, fc, mod_offset), 'omegaX', bx*g85) ...
 .set_pumping( R0 ) ...
 .solve().plot_solution();
